###   Implementation Constraints on natural climate solutions
###   Author: Kroeger, Timm; Erbaugh, James; Luo, Zhixian; Hilary Brumberg; Waverly Eichhorst; Margaret Hegwood; Anna LoPresti; Priya Shyamsundar; Peter W. Ellis; Lauren E. Oakes; Dow Martin; Pedro H. S. Brancalion; Mieke Bourne; Arundhati Jagadish; Kemen G. Austin; Andrew Kinzer; Marcos Sanjuán; Lisa McCullough; Marta Echavarria,

library(tidyverse)
library(ggplot2)
library(grid)
library(ggpattern)

nbase<-read.csv("NBaseData_v1.csv")


## Figure SI Heatmap by region

nbase <- nbase %>%
  mutate(NCS_Pathway = ifelse(NCS_Pathway == "regenerative agriculture (other than agroforestry)", "regenerative agriculturen\n (other than agroforestry)", NCS_Pathway)) %>%
  mutate(NCS_Pathway = ifelse(NCS_Pathway == "reduced woodfuel harvest in forests", "reduced woodfuel harvest\n in forests", NCS_Pathway)) %>%
  mutate(NCS_Pathway = ifelse(NCS_Pathway == "avoided coastal wetland conversion", "avoided coastal\n wetland conversion", NCS_Pathway)) %>%
  mutate(NCS_Pathway = ifelse(NCS_Pathway == "avoided grassland conversion", "avoided grassland\n conversion", NCS_Pathway))

subregion <- nbase %>% distinct(Subregion)

for (i in subregion$Subregion){
  nbase_heatmap <- nbase %>%
    filter(Subregion == i) %>%
    group_by(Constraint_category, NCS_Pathway) %>%
    summarise(frequency = n()) %>%
    ungroup()
  
  heatmapscaled <- nbase_heatmap %>%
    group_by(NCS_Pathway)%>%
    summarise(avg = mean(frequency), stDev = sd(frequency)) %>%
    full_join(nbase_heatmap,.)%>%
    mutate(., zscore= ((frequency-avg) / stDev)) %>%
    mutate(zscore_type = case_when(
      is.na(zscore) ~ "Undefined",
      TRUE ~ "Actual"
    ))
  
  allcombinations <- expand.grid(NCS_Pathway = unique(heatmapscaled$NCS_Pathway), Constraint_category = unique(heatmapscaled$Constraint_category))
  heatmapscaled <- merge(heatmapscaled, allcombinations, by = c("NCS_Pathway", "Constraint_category"), all.y = TRUE)
  heatmapscaled <- heatmapscaled %>%
    mutate(zscore_type = case_when(
      is.na(zscore_type) ~"NA",
      !is.na(zscore_type) ~ zscore_type
    ))
  
  scaled_heatmap<-ggplot(heatmapscaled, aes(x=NCS_Pathway, y=Constraint_category, fill=zscore))+
    geom_tile_pattern(aes(pattern = zscore_type), color = NA,
                      pattern_density = 0.1, pattern_spacing = 0.01, show.legend = TRUE) +
    scale_fill_gradient(name="Z-Scores", low = "white", high = "red", na.value = "grey")+
    scale_pattern_manual(values = c("Actual" = "none", "NA" = "none", "Undefined" = "stripe")) +
    coord_flip()+
    theme_classic()+
    ylab(NULL)+xlab(NULL)+
    theme(text = element_text(family = "Helvetica", size = 5),
          axis.text.x = element_text(angle = 45, hjust = 1),
          legend.text = element_text(""))+
    ggtitle(i) +
    guides(pattern = "none")
  
  scaled_heatmap
  
  file_name <- paste0("Fig S17 ",i, ".png")
  file_path <- paste0(file_name)
  ggsave(file_path, scaled_heatmap, width=88, height=88, units = "mm", dpi=300)
}